#include <cmath>
#include <iostream>
#include <vector>

using namespace std;

// 向量输出
void vec_out(vector<double> &u)
{
    cout << "(";
    for (int i = 0; i < u.size(); i++)
    {
        cout << u[i] << ",";
    }
    cout << "\b)";
}

// 矩阵输出
void mat_out(vector<vector<double>> &A)
{
    cout << "(" << endl;
    for (int i = 0; i < A.size(); i++)
    {
        vec_out(A[i]);
        cout << endl;
    }
    cout << ")";
}

// 矩阵转置
vector<vector<double>> trans_mat(vector<vector<double>> &A)
{
    vector<vector<double>> _A(A.front().size(), vector<double>(A.size()));
    for (int i = 0; i < A.size(); i++)
    {
        for (int j = 0; j < _A.size(); j++)
        {
            _A[i][j] = A[j][i];
        }
    }
    return _A;
}

// 获取矩阵的下三角部分
vector<vector<double>> lower_mat(vector<vector<double>> &A)
{
    vector<vector<double>> _A(A.size(), vector<double>(A.front().size()));
    for (int i = 0; i < A.size(); i++)
    {
        for (int j = 0; j < i + 1; j++)
        {
            _A[i][j] = A[i][j];
        }
    }
    return _A;
}

// 获取矩阵的上三角部分
vector<vector<double>> upper_mat(vector<vector<double>> &A)
{
    vector<vector<double>> _A(A.size(), vector<double>(A.front().size()));
    for (int i = 0; i < A.size(); i++)
    {
        for (int j = i; j < A.front().size(); j++)
        {
            _A[i][j] = A[i][j];
        }
    }
    return _A;
}

// 矩阵乘法
vector<vector<double>> mat_mult(vector<vector<double>> &A, vector<vector<double>> &B)
{
    vector<vector<double>> C(A.size(), vector<double>(B.front().size()));
    for (int i = 0; i < A.size(); i++)
    {
        for (int j = 0; j < B.front().size(); j++)
        {
            for (int k = 0; k < B.size(); k++)
            {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return C;
}

// 矩阵乘法
vector<double> mat_mult(vector<vector<double>> &A, vector<double> &b)
{
    vector<double> c(A.size());
    for (int i = 0; i < A.size(); i++)
    {
        for (int k = 0; k < b.size(); k++)
        {
            c[i] += A[i][k] * b[k];
        }
    }
    return c;
}

// 向量除法
void vec_div(vector<double> &u, double x)
{
    for (int i = 0; i < u.size(); i++)
    {
        u[i] /= x;
    }
}

// 向量内积
double vec_inner(vector<double> &u, vector<double> &v)
{
    double sum = 0;
    for (int i = 0; i < u.size(); i++)
    {
        sum += u[i] * v[i];
    }
    return sum;
}

// 求下三角阵的逆
vector<vector<double>> inv_lower(vector<vector<double>> L)
{
    vector<vector<double>> _L(L.size(), vector<double>(L.size()));
    for (int i = 0; i < L.size(); i++)
    {
        _L[i][i] = 1 / L[i][i];
        for (int j = 0; j < i + 1; j++)
        {
            L[i][j] /= L[i][i];
        }
    }
    for (int i=0;i<L.size();i++)
    {
        vector<vector<double>> tmp(L.size(), vector<double>(L.size()));
        for (int j=0;j<L.size();j++)
        {
            tmp[j][i] = -1 * L[j][i];
            tmp[j][j] = 1;
        }
        _L = mat_mult(tmp, _L);
    }
    return _L;
}

// Cholesky 分解
void cholesky(vector<vector<double>> &A)
{
    int n = A.size();
    double sum = 0;
    for (int i = 0; i < n; i++)
    {
        sum = 0;
        for (int p = 0; p < i; p++)
        {
            sum += A[i][p] * A[i][p];
        }
        A[i][i] = sqrt(A[i][i] - sum);
        for (int j = i + 1; j < n; j++)
        {
            sum = 0;
            for (int k = 0; k < i; k++)
            {
                sum += A[j][k] * A[i][k];
            }
            A[j][i] = (A[j][i] - sum) / A[i][i];
            A[i][j] = A[j][i];
        }
    }
}

// 求解上三角阵
vector<double> upper_tr(vector<vector<double>> &A, vector<double> &b)
{
    int m = A.front().size();
    int n = b.size();
    if (n > m)
    {
        n = m;
    }
    vector<double> x(n);
    double sum;
    for (int i = n - 1; i >= 0; i--)
    {
        sum = 0;
        for (int j = i + 1; j < n; j++)
        {
            sum += A[i][j] * x[j];
        }
        x[i] = (b[i] - sum) / A[i][i];
    }
    return x;
}

// 求解下三角阵
vector<double> lower_tr(vector<vector<double>> &A, vector<double> &b)
{
    int m = A.front().size();
    int n = b.size();
    if (n > m)
    {
        n = m;
    }
    vector<double> x(n);
    double sum;
    x[0] = b[0] / A[0][0];
    for (int i = 1; i < n; i++)
    {
        sum = 0;
        for (int j = 0; j < i; j++)
        {
            sum += A[i][j] * x[j];
        }
        x[i] = (b[i] - sum) / A[i][i];
    }
    return x;
}

// 返回指数幂
double pow(double x, int p)
{
    double res = 1;
    for (int i = 0; i < p; i++)
    {
        res *= x;
    }
    return res;
}

// 返回向量指数幂
vector<double> pow(vector<double> &x, int p)
{
    vector<double> y(x.size());
    for (int i = 0; i < x.size(); i++)
    {
        y[i] = pow(x[i], p);
    }
    return y;
}

// 求和内积
double inner(vector<double> u, vector<double> v)
{
    double sum = 0;
    for (int i = 0; i < u.size(); i++)
    {
        sum += u[i] * v[i];
    }
    return sum;
}

// 无穷范数
double inf_norm(vector<double> &x)
{
    double max = 0;
    for (int i = 0; i < x.size(); i++)
    {
        if (max < fabs(x[i]))
        {
            max = fabs(x[i]);
        }
    }
    return max;
}

// 最大位置
int max_pos(vector<double> &x)
{
    double max = 0;
    int index = 0;
    for (int i = 0; i < x.size(); i++)
    {
        if (max < fabs(x[i]))
        {
            max = fabs(x[i]);
            index = i;
        }
    }
    return index;
}

// householder 变换
double householder(vector<double> &x, vector<double> &v)
{
    double model = inf_norm(x);

    // protect x
    if (model == 0)
    {
        return 0;
    }

    vec_div(x, model);
    double sumOfx = vec_inner(x, x);
    double sumOfx1 = sumOfx - x[0] * x[0];
    v = x, v[0] = 0;
    double b;
    if (sumOfx1 == 0)
    {
        b = 0;
    }
    else
    {
        if (x[0] <= 0)
        {
            v[0] = x[0] - sqrt(sumOfx);
        }
        else
        {
            v[0] = -1 * sumOfx1 / (x[0] + sqrt(sumOfx));
        }
        b = 2 * v[0] * v[0] / (sumOfx1 + v[0] * v[0]);
        vec_div(v, v[0]);
    }
    return b;
}

// QR 分解
vector<double> QR(vector<vector<double>> &A)
{
    int row = A.size();
    int column = A.front().size();

    vector<double> d(column);
    for (int i = 0; i < column; i++)
    {
        vector<double> x(row - i);
        vector<double> v(row - i);
        for (int j = i; j < row; j++)
        {
            x[j - i] = A[j][i];
        }
        double b = householder(x, v);
        vector<double> w(column - i);
        // get w = b * AT * v
        for (int k = i; k < column; k++)
        {
            double sum = 0;
            for (int p = i; p < row; p++)
            {
                sum += A[p][k] * v[p - i];
            }
            // note: it is w[k-i], not w[k]
            w[k - i] = b * sum;
        }
        // get HA = A - v * wT
        for (int k = i; k < row; k++)
        {
            for (int p = i; p < column; p++)
            {
                if (p == i && k > i)
                {
                    A[k][p] = v[k - i];
                }
                else
                {
                    A[k][p] -= v[k - i] * w[p - i];
                }
            }
        }
        d[i] = b;
    }
    return d;
}

// 最小二乘法
vector<double> least_square(vector<vector<double>> &A, vector<double> &b)
{
    int m = A.size();
    int n = A.front().size();
    vector<double> d = QR(A);
    for (int i = 0; i < n; i++)
    {
        double sum = b[i];
        for (int j = i + 1; j < m; j++)
        {
            sum += A[j][i] * b[j];
        }
        sum *= d[i];
        b[i] -= sum;
        for (int j = i + 1; j < m; j++)
        {
            b[j] -= A[j][i] * sum;
        }
    }
    vector<double> x = upper_tr(A, b);
    return x;
}

// 幂法
double power(vector<vector<double>> &A)
{
    int n = A.size();
    vector<double> x(n, 1);
    double dp = 1, p, q = 0;
    vector<double> y(n);
    int index;
    while (dp > 1e-10)
    {
        y = mat_mult(A, x);
        index = max_pos(y);
        p = y[index];
        x = y;
        vec_div(x, p);
        dp = fabs(p - q);
        q = p;
    }
    return p;
}